package com.numericalmethod.algoquant.model.ralph2009.market;

import com.numericalmethod.suanshu.matrix.doubles.Matrix;
import com.numericalmethod.suanshu.matrix.doubles.matrixtype.dense.DenseMatrix;
import com.numericalmethod.suanshu.stats.stochasticprocess.multivariate.sde.MultivariateFt;
import com.numericalmethod.suanshu.stats.stochasticprocess.multivariate.sde.MultivariateSDE;
import com.numericalmethod.suanshu.stats.stochasticprocess.multivariate.sde.coefficients.DiffusionMatrix;
import com.numericalmethod.suanshu.stats.stochasticprocess.multivariate.sde.coefficients.DriftVector;
import com.numericalmethod.suanshu.vector.doubles.Vector;
import com.numericalmethod.suanshu.vector.doubles.dense.DenseVector;

/**
 * 
 * The model in "Ralph et al (2009) Momentum and Mean Reversion in Strategic Asset Allocation" 
 * consists of three stochastic processes: <br/>
 * 
 * <br/>
 * <li>dSt, equation 13 of Ralph 2009 paper: dSt/St = ( (mu0+mu1*Xt)*(1-phi) + phi * Mt ) dt + sigmaS1 dZt1 + sigmaS2 dZt2</li>
 * <li>dMt, equation 14 of Ralph 2009 paper: dMt = (1-phi)*(mu0 + mu1 * Xt - Mt) dt + sigmaS1 dZt1 + sigmaS2 dZt2</li>
 * <li>dXt, equation 15 of Ralph 2009 paper: dXt = -alpha * Xt dt + sigmaX1 dZt1 + sigmaX2 dZt2</li>
 * 
 * <br/>
 * We model them using the <code>MultivariateSDE</code> of SuanShu.
 * 
 * @author Paul/Stephen/Clement
 *
 */
public class Ralph2009MarketSDE extends MultivariateSDE {

	
	/**
	 * Create the system of SDEs mentioned in equation (13) (14) (15) in 
	 * "Ralph et al (2009) Momentum and Mean Reversion in Strategic Asset Allocation"
	 * 
	 */
	public Ralph2009MarketSDE(double mu0, double mu1, double phi, double alpha, 
				double sigmaS1, double sigmaS2, double sigmaX1, double sigmaX2) {
		
		// create a multivariate SDE with 2 independent driving Brownian motions
		super(new Ralph2009DriftVector(mu0, mu1, phi, alpha), 
				new Ralph2009DiffusionMatrix(sigmaS1, sigmaS2, sigmaX1, sigmaX2), 2);
		
	}
	

	
	
	/**
	 * Implementation of the drift term in 
	 * 
	 * <li>dSt, equation 13 of Ralph 2009 paper: dSt/St = ( (mu0+mu1*Xt)*(1-phi) + phi * Mt ) dt + sigmaS1 dZt1 + sigmaS2 dZt2</li>
     * <li>dMt, equation 14 of Ralph 2009 paper: dMt = (1-phi)*(mu0 + mu1 * Xt - Mt) dt + sigmaS1 dZt1  sigmaS2 dZt2</li>
     * <li>dXt, equation 15 of Ralph 2009 paper: dXt = -alpha * Xt dt + sigmaX1 dZt1 + sigmaX2 dZt2</li>
	 * 
	 */
	private static class Ralph2009DriftVector implements DriftVector {
		
		private double mu0;
		private double mu1;
		private double phi;
		private double alpha;

		public Ralph2009DriftVector(double mu0, double mu1, double phi, double alpha) {
			this.mu0=mu0;
			this.mu1=mu1;
			this.phi=phi;
			this.alpha=alpha;
		}
		
		@Override
		public Vector evaluate(MultivariateFt ft) {
			// interpret the filtration up to t.
			// expected to be 3 elements which are St, Mt, Xt in Ralph2009.
			Vector currentState=ft.Xt(); 
			if (currentState.size()!=3)
				throw new IllegalStateException("expecting filtration always have 3 state variables in Ralph2009DriftVector");
			
			double st=currentState.get(1);
			double mt=currentState.get(2);
			double xt=currentState.get(3);
			
			double sDrift = st * ((mu0 + mu1 * xt) * (1.0 - phi) + phi * mt);
			double mDrift = (1.0 - phi) * (mu0 + mu1 * xt - mt);
			double xDrift = -alpha * xt;
			
			return new DenseVector(sDrift, mDrift, xDrift);
		} 
	}
	
	/**
	 * Implementation of the diffusion term in 
	 * 
	 * <li>dSt, equation 13 of Ralph 2009 paper: dSt/St = ( (mu0+mu1*Xt)*(1-phi) + phi * Mt ) dt + sigmaS1 dZt1 + sigmaS2 dZt2</li>
     * <li>dMt, equation 14 of Ralph 2009 paper: dMt = (1-phi)*(mu0 + mu1 * Xt - Mt) dt + sigmaS1 dZt1  sigmaS2 dZt2</li>
     * <li>dXt, equation 15 of Ralph 2009 paper: dXt = -alpha * Xt dt + sigmaX1 dZt1 + sigmaX2 dZt2</li>
	 * 
	 */
	private static class Ralph2009DiffusionMatrix implements DiffusionMatrix {

		private double sigmaS1;
		private double sigmaS2;
		private double sigmaX1;
		private double sigmaX2;
		
		public Ralph2009DiffusionMatrix(double sigmaS1, double sigmaS2, double sigmaX1, double sigmaX2) {
			this.sigmaS1=sigmaS1;
			this.sigmaS2=sigmaS2;
			this.sigmaX1=sigmaX1;
			this.sigmaX2=sigmaX2;
		}
		
		/**
		 * Ralph2009 has 3 stochastic processes: St, Mt, Xt.
		 * 
		 * @return Hard-coded 3.
		 */
		@Override
		public int dimension() {
			// the 3 dimensions are namely St, Mt, Xt.
			return 3;
		}

		@Override
		public Matrix evaluate(MultivariateFt ft) {

			// interpret the filtration up to t.
			
			// expected to be 3 elements which are St, Mt, Xt in Ralph2009.
			Vector currentState=ft.Xt(); 
			if (currentState.size()!=3)
				throw new IllegalStateException("expecting filtration always have 3 state variables in Ralph2009DiffusionMatrix");
			
			double st=currentState.get(1);
			
			double[][] diffusionData = new double [3][];
			diffusionData[0]=new double [] { st * sigmaS1, st * sigmaS2 };
			diffusionData[1]=new double [] { sigmaS1, sigmaS2 };
			diffusionData[2]=new double [] { sigmaX1, sigmaX2 };
			
			return new DenseMatrix(diffusionData);

		}

		/**
		 * Ralph2009 has two independent brownian motion Zt1, Zt2
		 * 
		 * @return Hard-coded 2.
		 */
		@Override
		public int nB() {
			// 2 independent brownian motion namely Zt1, Zt2.
			return 2;
		}
		
	}
	
	
}
